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We perform a detailed analysis of the contact force network in a dense confined packing of pen- 
tagonal particles simulated by means of the contact dynamics method. The effect of particle shape 
is evidenced by comparing the data from pentagon packing and from a packing with identical char- 
acteristics except for the circular shape of the particles. A counterintuitive finding of this work is 
that, under steady shearing, the pentagon packing develops a lower structural anisotropy than the 
disk packing. We show that this weakness is compensated by a higher force anisotropy, leading to 
enhanced shear strength of the pentagon packing. We revisit "strong" and "weak" force networks in 
the pentagon packing, but our simulation data provide also evidence for a large class of "very weak" 
forces carried mainly by vertex-to-edge contacts. The strong force chains are mostly composed 
of edge-to-edge contacts with a marked zig-zag aspect and a decreasing exponential probability 
distribution as in a disk packing. 

PACS numbers: 



I. INTRODUCTION 



Among singular features of granular media, force trans- 
mission has received particular interest during the last 
decade. The contact forces in model granular media, as 
observed by experiments and numerical simulations, are 
highly inhomogeneous and their probability density func- 
tions (pdf's) are wide [J, i, i, i, B i, S i, H • The gran- 
ular texture is generically anisotropic in two respects: 
1) The contact normal directions are not random; 2) 
The force average as a function of contact normal di- 
rection is not uniform. The corresponding fabric and 
force anisotropies in shear are responsible for mechanical 
strength at the scale of the packing [l^, [111, [T^ , [Tsl . An- 
other interesting aspect, first analyzed in Ref. [lO| is the 
fact that the forces organize themselves in two distinct 
classes which contribute differently to fabric anisotropy, 
shear stress, and dissipation. In particular, the shear 
stress is fully transmitted via a "strong" contact net- 
work, materialized by force "chains" . The stability is 
ensured by the antagonist role of "weak" contacts which 
prop strong force chains fldi [T^ . 

The force transmission properties have been for the 
most part investigated in the case of granular media com- 
posed of isometric (circular or spheric) particles. How- 
ever, in various fields of science and engineering, the 
grains are seldom so "perfect". For example, elongated 
and platy shapes are encountered in biomaterials or phar- 
maceutical applications. Such shapes have unequal di- 
mensions and induce thus a degree of anisotropy in the 
bulk behavior in addition to fabric and force anisotropies 
[3) El) [13) Ell- On the other hand, granular geo- 
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materials are often composed of angular particles with 
plane faces as polyhedra. While rounded particles en- 
hance flowability, angular shape is susceptible to enhance 
the shear strength, a factor of vital importance to civil- 
engineering applications [H, [H, . The railway ballast 
is a well-known case where particle shape must be op- 
timized to avoid excessive differential settlement under 
vertical loading [2l|, [l^, . In such circumstances, the 
analysis of force transmission is a key to improve perfor- 
mance. 

In dealing with effects of particle shape, the issue is 
that a general quantitative description of particle mor- 
phology requires various shape parameters. For regular 
polygons in 2D, for instance, the only shape parameter 
is the number of sides (besides the diameter) whereas 
for irregular polygons more information is needed about 
the positions of the vertices in a reference system at- 
tached to the particle. In soil mechanics, angularity and 
roundedness are among basic parameters used to describe 
particle shapes d^]. As far as force transmission is con- 
cerned, at least two parameters seem to be most relevant: 
1) shape anisotropy (anisometry), which contributes to 
the anisotropy of stress transmission [ll|; 2) facettedness, 
which allows for extended (face to face, edge to face and 
edge to edge) contacts between particles leading possibly 
to the formation of columnar structures within a granular 
assembly. 

In this paper, we consider one of the simplest possi- 
ble shapes, namely regular pentagons. Among regular 
polygons, the pentagon has the lowest number of sides, 
corresponding to the least roundedness in this category, 
without the pathological space-filling properties of trian- 
gles and squares. We seek to isolate the effect of edge- 
to-edge contacts on force transmission by comparing the 
data with a packing of circular particles that, apart from 
the particle shape, is identical in all respects (prepara- 



2 



tion, friction coefRcients, particle size distribution) to the 
pentagon packing. Both packings are subjected to biaxial 
compression simulated by means of the contact dynam- 
ics method. The presence of edge-to-edge contacts affects 
both quantitatively and qualitatively the microstructure 
and the overall behavior during shear. These contacts do 
not transmit torques, but they are able to accommodate 
force lines that are usually unsustainable in packings of 
disks. 

This paper is organized as follows. We first present in 
Section |TT] the numerical procedures and a brief technical 
introduction to the detection and treatment of edge-to- 
edge contacts in the framework of the contact dynamics 
method. In Section lllli we compare stress-strain and 
volume-change characteristics. Then, In Sections IIVI and 
IVl we analyze the texture and force transmission fea- 
tures. In Section IVIl we focus on the pentagon pack- 
ing and we analyze the structure of force networks with 
vertex-to-edge and edge-to-edge contacts. The main re- 
sults are summarized and discussed in Section [VI 1 1 

II. NUMERICAL PROCEDURES 

The simulations were carried out by means of the con- 
tact dynamics (CD) method [H, HI]. The CD method 
is based on implicit time integration of the equations of 
motion and a nonsmooth formulation of mutual exclusion 
and dry friction between particles. This method requires 
no elastic repulsive potential and no smoothing of the 
Coulomb friction law for the determination of forces. For 
this reason, the simulations can be performed with large 
time steps compared to molecular dynamics simulations. 
We used the platform LMGC90 which is a multipurpose 
software developed in Montpellier, capable of modeling 
a collection of deformable or undeformable particles of 
various shapes [27| . 

A. Contact dynamics for polygons 

The particles are rigid polygons exerting normal and 
shear forces, /„ and ft, respectively, on each other. We 
attribute a positive sign to compressive normal forces. 
The relative normal velocity m„ between two particles in 
contact is counted positive when they move away from 
each other. Then, the condition of geometrical contact 
between two particles is expressed by the following mu- 
tually exclusive alternatives: 

fn ^ u„ = 
/„ = Un> 0. 

In the same way, the Coulomb friction law involves 
three mutually exclusive conditions: 

ft = -IJ-fn Ut > 

-^J'fn ^ ft ^ M/n Ut ^0 (2) 
ft = ^J■fn Ut <0 



where ut is the sliding velocity at the contact and n is 
the friction coefficient. The unknown variables are parti- 
cle velocities and contact forces. These are calculated at 
each time step by taking into account the conservation of 
momenta, the constraints expressed by ([T]) and ([2]), and 
the dissipation of kinetic energy during inelastic collisions 
between particles (ref). We use an iterative research al- 
gorithm based on a nonlinear Gauss-Seidel scheme. The 
uniqueness is not guaranteed for perfectly rigid particles 
in absolute terms. However, by initializing each step of 
calculation with the forces calculated in the preceding 
step, the set of admissible solutions shrinks to fiuctua- 
tions which are basically below the numerical solution. 
Let us note that in molecular dynamics simulations, this 
"force history" is encoded by construction in the particle 
positions. 

The research algorithm is applied to a set of potential 
contacts, identified or updated in each step. The contact 
detection between two bodies consists in looking for the 
overlaps of the portions of space they occupy. The treat- 
ment of the mechanical interaction requires additionally 
the identification of a common tangent plane (a line in 
2D). Of course, contact may take place through a larger 
contact zone than a single point. Several algorithms ex- 
ist for overlap determination between convex polygons 
[U, [23. In 2D simulations of the present paper, the de- 
tection of contact between two convex polygonal bodies 
was implemented through the so-called "shadow overlap 
method" devised by Moreau [ll], |23] , with reliability and 
robustness tested in several years of previo us appli cations 
to various states of granular materials [Tsl. l28l. [29j . 

In detection of contacts between two polygons, two sit- 
uations arise: 1) If a single corner is found crossing an 
edge of the partner polygon, the direction of this edge 
is viewed as the tangent direction. By orthogonally pro- 
jecting the intruding vertex onto the edge, one determines 
the penetration depth, while the nominal contact point 
is chosen at the center of this distance. Below, we will 
refer to this vertex-to-edge contact as "simple" contact. 
2) In case of double intrusion, the common tangent line is 
fixed from as a mean between the two overlapping edges 
and a segment of this line is identified as the contact 
segment. The impenetrability between two particles at 
such an edge-to-edge contact is ensured by applying the 
contact laws ([1]) and ^ to only two points of the con- 
tact segment (Fig. [T]). For this reason, we refer below 
to edge-to-edge contacts as "double" contacts. In prac- 
tice, two forces are calculated at each double contact, but 
only their resultant and application point are material. 
In this respect, the choice of the two points represent- 
ing a double contact does not affect the dynamics of the 
system. 



B. Numerical samples 

We generated two numerical samples. The first sample, 
denoted SI, is composed of 14400 regular pentagons of 
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FIG. 1: Representation of simple (vertex-to-edge) and double 
(edge-to-edge) contacts between two pentagons. 

three different diameters: 50% of diameter 2.5 cm, 34% of 
diameter 3.75 cm and 16% of diameter 5 cm. The second 
sample, denoted S2, is composed of 10000 discs with the 
same polydispersity. Both samples were prepared accord- 
ing to the same protocol. A dense packing was first con- 
structed following simple geometrical rules [13] and then 
compressed isotropically under a constant stress tro = lO'* 
Pa applied onto the right and top walls. The gravity was 
set to zero in order to avoid force gradients in the sam- 
ples. The coefficient of friction was set to 0.4 between 
grains and to with the walls. At equilibrium, both nu- 
merical samples were in isotropic stress state. The solid 
fraction was (j)Q = 0.80 for SI and (j)Q = 0.82 for S2. The 
aspect ratio was h/l k, 2, where h and I are the height 
and width of the sample, respectively. Figure [2] displays 
snapshots of the two packings at the end of isotropic com- 
paction. 

The isotropic samples were subjected to vertical com- 
pression by downward displacement of the top wall at 
a constant velocity of 1 cm/s for a constant confining 
stress (To acting on the lateral walls. The simulations 
were run up to a total cumulative vertical strain of 0.2 
with a time step of 5.10"'* s. The CPU time was 7.10"^^ 
s and 5.10"'* s per particle and per time step on a G5 
Apple computer. Since we are interested in quasistatic 
behavior, the shear rate should be such that the kinetic 
energy supplied by shearing is negligible compared to the 
static pressure. This can be formulated in terms of an 
"inertia parameter" / [3l[ defined by 

where e — y/y \s the strain rate, m is the total mass, 
and p is the average pressure. The quasistatic limit is 
characterized by the condition / <C 1. In our biaxial 
simulations, / was below 10"^. 

III. STRENGTH AND DILATANCY 

In this section, we compare the stress-strain and 
volume-change behavior between the packings of poly- 




FIG. 2: Snapshots of a portion of the samples S2 (a) and 
SI (b) composed of circular and pentagonal particles, respec- 
tively. 

gons (sample SI) and disks (sample S2). For the cal- 
culation of the stress tensor, we consider the "tensorial 
moment" of each particle i defined by [3, [s^l : 

Mlp^Y.r^r'p, (4) 

where is the a component of the force exerted on par- 
ticle i at the contact c, is the (3 component of the 
position vector of the same contact c, and the summa- 
tion is runs over all contacts c of neighboring particles 
with the particle i (noted briefly by c S i). It can be 
shown that the tensorial moment of a collection of rigid 
particles is the sum of the tensorial moments of individ- 
ual particles. The stress tensor cr for a packing of volume 
V is simply given by [TJ, [sll : 

where t'^ is the intercenter vector joining the centers of 
the two touching particles at the contact c. Remark that 
the first summation runs over all particles whereas the 
second summation involves all contacts in the volume 
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FIG. 3: Normalized shear stress g/p as a function of cumula- 
tive shear strain for the samples SI and S2. 



V, with each contact appearing only once. We extract 
the mean stress p — (cti + cr2)/2, and the stress devia- 
tor q = ((Ti — (T2)/2, where cti and (T2 are the principal 
stresses. The major principal direction during vertical 
compression is vertical. 

The strain parameters are the cumulative vertical, hor- 
izontal and volumetric strains ei, £2 and £p, respectively. 
By definition, we have 



£1 



dh _ / Ah 
h \ ho 



(6) 



where ho is the initial height and Ah — ho — his the total 
downward displacement, and 



dV 



= In 1 



Ay 

Vo 



(7) 



where Vo is the initial volume and AV = V — Vo is the 
cumulative volume change. 

Figure [3] shows the normalized shear stress q/p for the 
samples SI and S2 as a function of shear strain Sq = Si~ 
£2- For S2, we observe a classical behavior characterized 
by a hardening behavior followed by (slight) softening 
and a stress plateau corresponding to the residual state 
of soil mechanics For SI, we observe no marked 

stress peak. The residual stress is higher for polygons 
(~ 0.35) than for disks (~ 0.28). This means that the 
polygon packing has a higher angle of internal friction (p 
defined by 



sm If = —. 
P 



(8) 



Figure [4] displays the cumulative volumetric strain £p 
for polygons and disks as a function of Eg. Both sam- 
ples dilate and tend to isochoric deformation at large 
strains. It is remarkable that the polygon packing SI 
initially dilates less than the disk packing S2. This be- 
havior is reversed at larger strains with a crossover occur- 
ring after the peak state. Notice that the solid fraction 
is initially lower in SI (0.80) than in S2 (0.82). This is 
because it is more difficult to obtain a compact pack- 
ing with polygonal shapes by isotropic compression as 
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FIG. 4: Cumulative volumetric strain Sp as a function of cu- 
mulative shear strain for the samples SI and S2. 
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FIG. 5: Stress-dilatancy relation between dilatancy angle tp 
and internal angle of friction ifi for the samples SI and S2. 



a result of enhanced steric effects compared to disks. In 
other words, angular particles can form larger pores com- 
pared to rounded particles. The volumetric deformation 
can also be expressed in terms of the so-called " dilatancy 
angle" ip defined by [s^ 



sin ip 



(9) 



The cumulative angle of dilatancy, i.e. during shear up to 
the residual state, is only slightly higher for the polygon 
packing than the disk packing. 

The plot of i/j as a function of ip, i.e. the so-called 
stress-dilatancy diagram, is shown in Fig. [5] for polygons 
and disks [33,] . Remarkably, both plots are parallel to the 
line if = ip with an offset ipo'- 



ip2^ ipO+Tp. 



(10) 



The offset (po is the friction angle at zero dilatancy. We 
have (po — 0.12 for disks and ipo — 0.3 for polygons. 
This observation is in agreement with the arguments of 
Taylor [33, 34] based on energy balance and recently re- 
visited also in the case of cohesive granular media [35j . 
The higher level of ip for the polygon packing reflects 
the organization of the microstructure and the features 
of force transmission for each particle shape. This point 
is considered in more detail in the following section. 
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IV. GRANULAR TEXTURE 

The granular texture, i.e. the organization of the parti- 
cles and their contacts in space, is basically controlled by 
steric exclusions between the particles and force balance 
conditions [s^. The texture can be described in terms 
of various statistical descriptors pertaining to the force- 
bearing network of particles. At the lowest order, the 
compactness of the structure can be described in terms 
of both the solid fraction p and the coordination number 
z. The connectivity of the network can further be char- 
acterized by the fraction P(c) of particles having exactly 
c contact neighbors. These are scalar parameters or func- 
tions. At higher orders, the anisotropy of the texture is 
described by different "fabric tensors" . We consider here 
these geometrical descriptors in order to identify the sig- 
nature of particle shape. 



A. Connectivity 

The connectivity of the particles by force-bearing con- 
tacts is described at the lowest order by the average num- 
ber z of contact neighbors per particle. The particles 
with no force-bearing contact are thus removed from the 
statistics. Note also that each double (edge-to-edge) con- 
tact for the polygons is counted once although double 
contacts are treated as two point contacts belonging to 
the contact segment (see section . Fig. |6K displays 
the evolution of z for the pentagon packing (SI) and the 
disk packing (S2) as a function of Sq. The coordination 
number evolves to a steady-state value in both samples 
that is higher for S2 (~ 3.85) than for SI (~ 3.75). The 
difference is, however, much less important than in the 
initial configuration (~ 3.95 for S2 compared to ~ 3.20 
for SI) prepared by means of isotropic compaction. 

It is also interesting to compare the two samples in 
terms of "contact lifetimes" . Let us consider a reference 
configuration, e.g. the initial state of each sample. We 
follow the history of each contact listed in this state. In 
particular, we define 7 as the fraction of persistent con- 
tacts of the initial list. During deformation, 7 declines 
from 1 to as an increasing number of initial contacts 
are lost due to particle rearrangements. Fig. [6|d shows 
gamma as a function of Eq for SI and S2. We see that, 
following a rapid initial falloff, 7 decreases slowly in both 
samples but the rate of contact loss is globally higher for 
polygons than disks. We remark that even at Sq — 0.4, 
the contact list is renewed by only 50%. 

The connectivity P{c) of the particles is plotted in Fig. 
[7] for SI and S2 at £q = 0.3. Interestingly, the two plots 
are nearly identical with a peak for c — A. In both sam- 
ples, the fraction of particles with 5 contacts is larger 
than that with 3 contacts. This shows that the connec- 
tivity does not reflect the difference in texture between 
the two packings although a qualitative difference exists 
as we shall see below by considering fabric anisotropy and 
force transmission. 



4.0 

3.8 
3.6 
3.4 
3.2 











(a) 




— SI 
-- S2 












0. 



FIG. 6: The coordination number z (a) and the fraction 7 
of persistent contacts (b) as a function of cumulative shear 
strain Sq for the samples SI and S2. 
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FIG. 7: Connectivity diagram for the samples SI and S2 ex- 
pressing the fraction P{c) of particles with exactly c contacts 
in the residual state. 



B. Fabric anisotropy 

The shear strength of dry granular materials is gener- 
ally attributed to the buildup of an anisotropic structure 
during shear due to friction between the particles and 
as a result of steric effects depending on particle shape 
i37j 38, 39J. Several methods have been used to quantify 
the fabric (structural) anisotropy of granular materials 
(T3I [40l . . A common approach is to consider the prob- 
ability distribution P{n) of the contact normals n which 
are generically nonuniform. In two dimensions, the unit 
vector n is described by a single angle 9, the orientation 
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FIG. 8: Polar representation of the probability density func- 
tion Pe of the contact normal directions 9 for the samples SI 
and S2 in the residual state. 

of the contact normal. The probability density function 
Pg{0) of contact normals provides a detailed statistical 
information about the fabric. It is 7r-periodic in the ab- 
sence of an intrinsic polarity for n. 

Most structural information is generally condensed in 
the second moment of Pe, called fabric tensor [40j : 

Fo.p = l r no.{e)n0{e)Pg{e)de ^ j-Y,nl,n% (11) 

^ CSV 

where a and /? design the components in a reference 
frame and Nc is the total number of contacts in the con- 
trol volume V . By definition, tr{F) = 1. The anisotropy 
of the contact network is given the difference between 
the principal values Fi and F2. We define the fabric 
anisotropy a by 

a^2iFi-F2). (12) 

For fix coordinates, with the x-axis pointing along 0' , we 
define also a "signed anisotropy" a' by 

a' ^2{Fi- F2)cos2{0F ~0'), (13) 

where 9p is the major principal direction of the fab- 
ric tensor. For 6' = Op, we have a' = a. The signed 
anisotropy corresponds to the second term of the Fourier 
expansion of Pe{0) and it is useful whenever the direction 
of anisotropy is not constant. 

Figure [8] displays a polar representation of Pq (9) for 
the samples SI and S2 at = 0.3. We observe a 
nearly isotropic distribution for the pentagon packing in 
spite of shearing whereas the disk packing is markedly 
anisotropic. This is a surprising observation in view of 
the higher shear strength of the pentagon packing (Fig. 
|3]). It is also counterintuitive as one expects that double 
contacts should allow a polygon packing to build more 
easily an anisotropic structure. 

The evolution of a' is shown in Fig. [5] as a function 
of Eq for SI and S2. The privileged direction of the con- 
tacts, corresponding to Op, is vertical in both packings. 
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FIG. 9: Evolution of the anisotropy a' with cumulative shear 
strain Sq for the samples SI and S2. 

In both cases, a' increases from (as a result of the initial 
isotropic compression) to its largest value in the residual 
state. The anisotropy stays quite weak in the pentagon 
packing whereas the disk packing is marked by a much 
larger anisotropy, increasing to ~ 0.3 and then relaxing 
to a slightly lower value in the residual state. As we shall 
see below, the low anisotropy of the pentagon packing re- 
sults from a particular organization of the force network 
in correlation with the orientations of simple and double 
contacts in the packing (section |VI| . We will also show 
that the large shear strength of the pentagon packing is a 
consequence of a strong force anisotropy in this packing 
(see next section). 

V. FORCE TRANSMISSION 

In this section, we analyze the anisotropy and inhomo- 
geneity of force networks in the packings of pentagons 
and disks. This leads us to consider the contributions of 
force and texture anisotropics to average shear stresses. 



A. Force anisotropy 

The angular distribution of contact forces in a granular 
packing can be represented by the average force {f){n) 
as a function of the contact normal direction n. We dis- 
tinguish the average normal force (/„) from the average 
tangential force {ft) formally defined by [l^ 

{ftm - T^) T ft, ^ ^ 

where and are the normal and tangential forces, 
respectively, acting at the contact c (according to a 
sign convention attributing positive values to the nor- 
mal forces), S{9) is the set of contacts with direction 
9 £ [8 - Ae/2,9 + A9/2] for angle increments A9, and 
Nc{9) is the number of contacts in S{9). 



FIG. 10: Polar representation of tfie angle- averaged normal ' ^ ^ ^ 

(a) and tangential (b) forces {fn){0) and {ft){d) for the sam- % 0.1 0.2 0.3 0.4 

pies SI and S2 in the residual state. £ 



By definition, the two functions (/„) and (ft) are tt- 
periodic. After sufficiently long monotonous shearing, 
these functions can be approximated by their Fourier ex- 
pansions truncated beyond the second term [3, ■ 

j {fn){9) - (/){l + a„cos2(0-0„)} 
I {fm - {f)at8m2{e-0t), 

where (/} is the average force, a„ and at represent the 
anisotropics of the normal and tangential forces, respec- 
tively, and 6n and 9t are their privileged directions. 

In Fig. \M the functions (/„)(6') and {ft) (9) are dis- 
played in polar coordinates at Sq = 0.3. The pentagon 
and disk packings show pronounced force anisotropy with 
a stronger anisotropy in the case of pentagons both for 
normal and tangential forces. These plots can be fitted 
by harmonic functions [Eq. (jl5[) ] in order to estimate the 
force anisotropics a„ and at- However, it is more conve- 
nient to estimate the anisotropics through the following 
"force tensors" : 

f = J{fn){0)n^npde, 

I 'i (16) 

= lifmn^npde. 
K. 

It is easy to see that tr{H^"^) = tr{H^*'>) = (/), and by 
identification with ([T5)) we have 











rr(t) 









(17) 



where the subscripts 1 and 2 refer to the principal values 
of the tensors. 



FIG. 11: Evolution of force anisotropics a„ (a) and at (b) as 
a function of cumulative shear strain Sq in samples SI and S2. 



Figure [Til shows the evolution of a„ and at with Sq in 
samples SI and S2. We see that, in contrast to fabric 
anisotropics (Fig. [9]), the force anisotropics in pentagon 
packing remain always above those in the disk packing. 
This means that the aptitude of the pentagon packing to 
develop large force anisotropy and strong force chains is 
not solely dependent on the global fabric anisotropy of 
the system. In section IVI[ we will show that the force 
anisotropy of the pentagon packing stems from the high 
anisotropy of the sub-network of double contacts and 
strong activation of friction forces. Indeed, due to the 
geometry of the pentagons, i.e. the absence of parallel 
sides, the strong force chains are mostly of zig-zag shape, 
as observed in Fig. [T3b . and the stability of such struc- 
tures requires strong activation of tangential forces. This 
explains, in turn, the large value of at for pentagons, very 
close to a„, whereas in the disk packing at is nearly half 
of a„. 

The anisotropics a, a,i and at are interesting descrip- 
tors of granular microstructure and force transmission as 
they underlie the shear stress. Indeed, it can be shown 
that the general expression of the stress tensor Eq. ^ 
leads to the following simple relation [l^, [3§| : 

- ~l-{a + an + at), (18) 
P 2 

where the cross products between the anisotropics have 
been neglected and it has been assumed that the stress 
tensor is coaxial with the fabric tensor Eq. ([TTjl and the 
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force tensors Eq. Fig. [T^ shows that Eq. HH]) holds 

quite well both for pentagons and disks. 

A remarkable consequence of Eq. ([T5)) is to reveal the 
distinct origins of shear stress in pentagon and disk pack- 
ings. The fabric anisotropy provides a major contribution 
to shear stress in the disk packing (Fig. [5]) whereas the 
force anisotropics are more important for shear stress in 
the pentagon packing (Fig. [Tl]). In this way, in spite of 
the weak fabric anisotropy a, the larger force anisotropics 
a„ and at allow the pentagon packing to reach higher lev- 
els oi q/p compared to the disk packing. 



B. Force distributions 




FIG. 13: (color online) Snapshots of normal forces in samples 
S2 (a) and SI (b). Line thickness is proportional to the normal 
force. 



The strong inhomogeneity of contact forces is a well- 
known feature of granular media. It has been investigated 
mostly for spherical or cylindrical particles both by ex- 
periments and numerical simulations P, H, 0, H, d, 0, B @ ■ 
The probability density function (pdf ) of normal forces is 
characterized by two features which seem to be specific 
to granular media: 1) The pdf is roughly a decreasing 
exponential function for forces above the mean, 2) In the 
range of weak forces below the mean, the pdf does not 
decline to zero with the force. The relative scatter of data 
reported by different authors for weak forces shows the 
sensitivity of the pdf in this range to the details of the mi- 
crostructure. But, the common observation that there is 
a large number of contacts transmitting very weak forces, 
is a straightforward signature of the arching effect. From 
this point of view, one expects that angular particle shape 
will influence mainly the distribution of weak forces by 
enhancing the arching effect. 

Figure [T2] displays maps of normal forces in a portion 
of each of the samples SI and S2 at a large cumulative 
strain. We observe the strong anisotropy of normal forces 
in the pentagon packing compared to the disk packing (as 
discussed in section |V]) as well as the zig-zag form of the 
strong force chains. The normal force pdf's are shown in 



Fig. [13] in log-linear and log-log scales at large strains. 
The forces are normalized by the mean normal force (/„) 
in each sample. In both samples, the number of strong 
forces (above the mean (/„)) falls off exponentially: 



/„ oc e-"!-'^"/^/") in SI, 
/„ oc e-"2-^"/</"> inS2, 



(19) 



with ai ~ 0.74 and a2 — 1.4. The smaller value of 
ai means that the distribution is wider for pentagons 
compared to disks. The distribution is nearly uniform in 
the whole range of weak forces (/„ < (/n)) in S2. In the 
pentagon packing SI, we observe a uniform distribution 
only in the range 0.1(/„) < fn < ifn)- Nearly 30% of 
forces are in this range. The number of "very weak" 
forces in SI in the range /„ < 0.1(/„) increases faster 
than a power law as /„ tends to zero. A fraction ~ 
30% of contacts belong to this range. The presence of 
numerous "very weak" forces in the pentagon packing is 
a clear signature of enhanced arching effect that can be 
characterized, as we shall see below, by the respective 
roles of simple and double contacts with respect to force 
transmission. 
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FIG. 14: Probability density functions of normal forces in 
samples SI and S2 in log-linear (a) and log-log scales (b). 



C. Bimodal character of stress transmission 

The genuine organization of contact forces in granu- 
lar media, involving strong force chains propped by weak 
forces, was first analyzed by Radjai et al. by means of 
contact dynamics simulations for packings of circular and 
spherical particles [Ifl] . This analysis proceeds by consid- 
ering the subset of contacts which carry a force below a 
cutoff force ^ normalized by the mean force. This sub- 
set is referred to as the "^-network" . The variation of 
a quantity evaluated for the "^-network" as ^ is varied 
from to the maximal force in the system, provides its 
correlation with the contact force. Here, we apply this 
same approach to SI and S2 samples for the stress ratio 
i{C)Ip-i defined as stress deviator q(£_) (normalized by the 
total pressure p of the sample) in the ^-network, and for 
a(^), defined as the fabric anisotropy in the ^-network. 

The plot of q{C}/p is shown in Fig. [T5]for SI and S2 in 
the residual state. In both samples, the stress deviator is 
nearly zero for ^ < 1, i.e. for the normal forces below the 
average force. This means that the shear stress is almost 
totally sustained by the "strong" contact network ^ > 1 
for the pentagon packing as well as for the disk packing. 
Fig. [TBI shows the fabric anisotropy a'(^) as a function 
of ^ in the samples SI and S2. By definition, a positive 
value of a' corresponds to the principal stress direction 
whereas a negative value corresponds to the orthogonal 
direction. We see that the direction of anisotropy is or- 
thogonal to the principal stress direction (a' < 0) for 




FIG. 15: Partial shear stress q{C)/p as a function of force 
cutoff 5 (normalized by the mean force) for the samples SI 
and S2 in the residual state. 




FIG. 16: Partial fabric anisotropy a'(^) as a function of force 
cutoff ^ (normalized by the mean force) in the samples SI and 
S2. 



weak forces (small f). This "orthogonal" anisotropy of 
the weak forces is more important in the pentagon pack- 
ing compared to the disk packing, and, as shown in the 
inset to Fig. 1161 it is mainly due to "very weak" forces. 
When ^ is increased beyond (/„), a' becomes less neg- 
ative and finally changes sign, showing that the strong 
contacts are preferentially parallel to the principal axis. 
These strong contacts are less than 40% of all contacts, 
but their positive contribution to a' overcompensates the 
negative contribution weak contacts. For large ^, the par- 
tial anisotropy approaches the fabric anisotropy of the 
whole system. 

These data demonstrate the bimodal character of 
stress transmission also in the pentagon packing in spite 
of a very different particle geometry. The mean force 
plays a particular role in differentiating strong contacts 
from weak contacts. However, the force pdf's (Fig. 
rn]) and the anisotropy of weak forces (Fig. [TB)) pro- 
vide also evidence for the existence of a class of very 
weak forces, corresponding approximately to the range 
/„ < 0. !(/„), within the weak network. This class is 
strongly anisotropic with a privileged direction which is 
orthogonal to the major principal stress direction, and 
the corresponding force pdf diverges as the force tends 
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FIG. 17: (color online) Tricolor map of the contact net- 
work composed of very weak (blue) , intermediate (green) and 
strong (red) contacts in the pentagon packing. 
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FIG. 18: Normalized shear stress q/p for simple (s) and double 
(d) contacts, as well as for all contacts (s+d), as a function 
of cumulative shear strain Sq in the pentagon packing. 

to zero. Fig. [T7] displays a tricolor map of the contact 
network representing very weak, intermediate (0.1(/„) < 
fn < (/«)) and strong contacts in the pentagon packing. 
Large cells of strong contacts are composed of zig-zag 
chains. The anisotropy of strong contacts is reflected in 
the elongated shape of these cells along the major prin- 
cipal stress direction. Both intermediate and very weak 
forces prop these cells. 



VI. SIMPLE VERSUS DOUBLE CONTACTS 

In this section, we focus on the organization of simple 
and double contacts in the pentagon packing. The double 
contacts, i.e. the side-sharing polygons, are generally as- 
sumed to be at the source of the higher strength of poly- 
gon packings. For the texture, we would like also to inves- 
tigate the proportions of simple and double contacts and 
their respective contributions to the overall anisotropy of 
the pentagon packing. It is also important to identify the 
role of double contacts in force transmission. 
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FIG. 19: Proportions ks and kd of simple and double contacts, 
and the corresponding relative force averages fs and fd, as a 
function of cumulative shear strain Eq. 

The general expression [Eq. [5] of the stress tensor cr 
allows us to perform a unique additive decomposition of 
the stress into two parts: 

cr = <T, + o-d, (20) 

where (Tg is obtained from the expression ([5]) by restrict- 
ing the summation to simple contacts, and cr^ is the com- 
plementary tensor involving only double contacts. The 
respective stress deviators qg and qd normalized by the 
mean stress p are shown in Fig. [18] as a function of 
strain £q. The strength qd/p of double contacts varies 
from two to three times that of simple contacts during 
shear deformation of the pentagon packing. The propor- 
tions ks and kd of simple and double contacts are shown 
in Fig. [in] as a function of Eq. The same figure dis- 
plays the relative force averages fs = ks{fn)d/ (fn) and 
fd = kd{fn)s/{fn), where (/„)« and {fn)d are the mean 
normal forces of simple and double contacts, respectively. 
We see that kd increases with strain but remains below 
ks- On the other hand, initially we have fd = fs = 0.5, 
reflecting the isotropic state of the packing prepared by 
isotropic compaction. However, fd increases with shear 
up to fd — 1.5/s in the residual state. This means that 
the larger shear stress carried by double contacts in the 
residual state is due to the larger mean normal force of 
double contacts despite their smaller proportion in the 
packing. 

The growth of the number of double contacts shown 
in Fig. [19] represents the gradual consolidation of the 
sample. In Fig. [2D]we plot the cumulative proportions 
Ajs^d and Ajd^s of simple contacts turning to double 
and vice versa, respectively. Although transformation 
between the two contact types occurs at each step in both 
directions s — > d and d s, the consolidation involves 
on average a net fraction of simple contacts transforming 
into double contacts. 

The connectivity of the pentagon packing by simple 
and double contacts can be represented by the proportion 
P{ms,md) of particles with exactly TOs simple contacts 
and rud double contacts. Fig. [^l] shows a grey level map 
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FIG. 20; Cumulative proportion A7 of simple contacts turn- 
ing to double (s d) and vice versa (d — > s). 
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FIG. 21: Grey level map of the connectivity function 
P{ms,md) of the pentagon packing in the residual state. 

of this function for the pentagon packing in the residual 
state. The row rud = corresponds to particles with only 
simple contacts (nearly 2% of the total number of parti- 
cles) whereas the columnm^ = represents the particles 
with only double contacts (nearly 6%). On average, a 
particle has more simple contacts than double contacts 
but the maximum occurs at rus — nid — 2. 

We now consider the fabric tensor decomposed in a 
similar way as the stress tensor [Eq. ((20|) ] into two partial 
tensors: 

F = Fs+Fd, (21) 

where Fg and Fd are defined as F in Eq. pT|) by simply 
restricting the summation to simple and double contacts, 
respectively, and by dividing the sum by the total num- 
ber Nc of contacts. The respective anisotropics a'^ and 

of simple and double contacts are displayed in Fig. 

as a function of Eg. The interesting observation here 
is that the simple contacts have a negative anisotropy 
which, according to Eq. , means that simple contacts 
are mostly oriented perpendicular to the major principal 
fabric direction Op. In other words, most simple con- 
tacts belong to the weak network. In contrast, the double 
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FIG. 22: The anisotropy a of simple (s) and double (d) con- 
tacts as a function of cumulative shear strain Sq in the pen- 
tagon packing. 



contacts have an increasing positive anisotropy which is 
larger than the mean anisotropy a of the sample. This 
is consistent with the fact that the double contacts take 
over larger forces and they contribute more to the shear 
stress than simple contacts. 

The normal force pdf 's for simple and double contacts 
are shown in Fig. [221 Both contact types are involved in 
weak and strong networks and the pdf's have the same 
functional form. But the contribution of simple contacts 
is more important in the range of very weak forces. Once 
again, as for anisotropy, the very weak contacts appear to 
be related to the particular geometry of the pentagons. 
At large strains, about 32% of all contacts belong to 
the very weak force network with 25% simple contacts 
against 7% double contacts. A snapshot of the normal 
force network is shown in Fig. [24] where the line widths 
are proportional to the line width with different colors 
(or grey levels) for simple and double contacts. The re- 
markable feature of this map is the network of very strong 
zigzag force chains composed mostly of double contacts 
and occasionally mediated by simple contacts. 

The proportions fcf and of strong (S) and weak 
(W) simple (s) contacts, respectively, as well as the pro- 
portions and of strong and weak double (d) con- 
tacts are plotted in Fig. [25|as a function of £q. We see 
that in the strong network (/„ > (/)) the proportion 
of double contacts is nearly the same as the proportion 
fcf of simple contacts in the initial (isotropic) state, but 
during shear fcf declines down to kf ~ 0.5fc^ in the resid- 
ual state, in agreement with the impression left by Fig. 
[231 We have an inverse situation for the weak network 
composed of two times more simple contacts than double 
contacts, i.e. kY — ^k^ in the residual state. It is also 
interesting to remark that the fraction of weak contacts, 
i.e. + kY — 0.58 in the residual state is very close to 
that (0.62) in the case of the disk packing. 



12 




10' 



T3 



10 



10' 



10" 



10 



10 



10 



FIG. 23: Probability density function of normal forces for 
simple (s) and double (d) contacts in log-linear (a) and log- 
log scales (b). 




FIG. 24: (color online) Color map of the normal force network 
in the residual state with simple contacts (s) in blue and dou- 
ble contacts (d) in red. Line thickness is proportional to the 
normal force. 



VII. CONCLUSION 

The objective of this paper was to isolate the effect of 
particle shape on force transmission in granular media by 
means a detailed comparison between two similar pack- 
ings with different particle shapes: pentagons vs. disks. 
We observed enhanced shear strength and force inhomo- 
geneity in the pentagon packing. But, unexpectedly, the 




FIG. 25: Proportions fcf and of strong (S) and weak (W) 
simple (s) contacts, respectively, as well as the proportions fcf 
and kY of strong and weak double (d) contacts as a function 
of cumulative shear strain Sq in the pentagon packing. 



pentagon packing was found to develop a lower struc- 
tural (fabric) anisotropy compared to the disk packing 
under shear. This low fabric anisotropy, however, does 
not prevent the pentagon packing from building up a 
strong force anisotropy that underlies its enhanced shear 
strength compared to the disk packing. 

This finding is interesting as it shows unambiguously 
that the force anisotropy in a granular material has two 
distinct sources: (1) Fabric anisotropy, with a maximum 
value depending on particle shapes; (2) Particle shapes. 
The first mechanism is crucial for the disk packing so 
that the force anisotropy, and the shear stress as a re- 
sult, vanishes in an isotropic disk packing (e.g. when the 
friction coefficient is set to zero) . The second mechanism 
may be the predominant source of strength for "facetted" 
particles that can give rise edge-to-edge (in 2D) or face- 
to-face (in 3D) contacts allowing for strong force localiza- 
tion along such contacts in the packing. Since the fabric 
anisotropy is low in a pentagon packing, the role of force 
anisotropy and thus the local equilibrium structures or 
arching are important with respect to its overall strength 
properties. The pentagons analyzed in this work provide 
thus the first counter-example of a system where the role 
of fabric anisotropy in shear strength is marginal. 

Another shape- related effect was the observation of zig- 
zag force chains mostly composed of edge-to-edge con- 
tacts in steady shearing. The vertex-to-edge contacts 
belong thus mainly to the weak force network or a class 
of "very weak" forces that can be considered as a signa- 
ture of enhanced arching or screening effect of forces in 
the presence of edge-to-edge chains. These "very weak" 
forces can also be observed, though to a lesser extent, 
in a disk packing with high coefficients of friction Q 
or on experimental pdf 's of normal forces acting on the 
walls of a container [5|]- Let us recall that a "very weak 
phase" was also evidenced by considering the correlation 
between friction mobilization and the anisotropy of gran- 
ular texture in a disk packing at the stability limit |14l | . 

By focusing on pentagon packings, we were able 
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to demonstrate the nontrivial phenomenology resulting 
from the specific shape of particles as compared to a disk 
packing. Although general features of force transmis- 
sion (pdf's, bimodal character, etc) seem to be robust, 
the details of force transmission (relative importance of 
force and fabric; anisotropy, the role of edge-to-edge con- 
tacts, etc) seem to be strongly shape-dependent. Cur- 



rently, wc work to elucidate this issue for regular poly- 
gons (hexagons and higher number of sides) as well as 
polyhedral particles in three dimensions. 
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SNCF, and the Region Languedoc-Roussillon of France. 
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